Fermi-Pasta-Ulam (3 lattice: Peierls equation and anomalous heat conductivity 
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The Peierls equation is considered for the Fermi-Pasta-Ulam lattice. Explicit form of the 
linearized collision operator is obtained. Using this form the decay rate of the normal mode energy 
as a function of wave vector k is estimated to be proportional to fc 5//3 . This leads to the t~ 3 ^ 5 long 
time behavior of the current correlation function, and, therefore, to the divergent coefficient of heat 
conductivity. These results are in good agreement with the results of recent computer simulations. 
Compared to the results obtained through the mode coupling theory our estimations give the same 
k dependence of the decay rate but a different temperature dependence. Using our estimations we 
argue that adding a harmonic on-site potential to the Fermi-Pasta-Ulam f3 lattice may lead to finite 
heat conductivity in this model. 

PACS numbers: 63.10.+a, 05.60.-k, 44.10.-M, 66.70.+f 
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I. INTRODUCTION 

The Peierls equation has played an important role 
in understanding properties of solids since its original 
derivation by Peierls 0, Q . It was successfully used for 
qualitative explanation of heat conduction in dielectrics 
and for prediction of such phenomena as second sound 
and Poiseuille flow Q . In spite of these successes quanti- 
tative predictions are hard to make due to the enormous 
complexity of the equation even for solids with simple 
dispersion laws. It is well known that approximating the 
solid by isotropic continuum leads to divergent heat con- 
ductivity even in three dimensions if only three phonon 
collisions are considered 0, • This divergence can be 
eliminated if more careful analysis of the dispersion re- 
lations is performed In this paper we would like to 
consider the collision operator of the linearized Peierls 
equation for a simple one dimensional model: linear chain 
with quartic interaction also known as Fermi-Pasta-Ulam 
(FPU) P lattice. 

One dimensional lattices has drawn considerable atten- 
tion since the original work of Fermi, Pasta, and Ulam 
@ and originated a vast area of research Q. It is now 
well established by computer simulations that the heat 
conductivity in FPU lattices diverges when the size of the 
lattice goes to infinity |a, |£| |lfj, llll Il2| . The simulations 
give a power law dependence of the heat conductivity on 
the number of particles N as approximately N 2 / 5 . This 
form of N dependence is related to the i -3 / 5 long time 
behavior of the current correlation function. Theoretical 
work in this area 0, IT(i| was based on the application of 
the mode coupling theory 0, Q . Application of Peierls 
equation to the analysis of heat conduction in FPU lat- 
tices was limited to qualitative estimations pi 0, 0] . It 
is therefore of interest to check if more careful analysis of 



the Peierls equation can explain some of the anomalous 
properties of the FPU chains. 

To this end we will consider the explicit form of the 
Peierls equation for the FPU lattice in Sec. II. In Sec. Ill 
we apply this equation to estimate the long time behavior 
of the current correlation function as well as the wave 
vector dependence of the decay rate for mode energies. 
Concluding remarks are given in Sec. IV. 



II. THE PEIERLS EQUATION 



The Hamiltonian for the FPU j3 lattice is 
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Here u r is the displacement of the particle at site r, p r is 
the momentum conjugate to u r and m is the mass of the 
particle. We also use A as a coupling constant. Cyclic 
boundary conditions are imposed, i.e., u r+ N = u r , where 
N is the number of particles. We can introduce action 
variables Jk and angle variables ctk related to u r and p r 
through 
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Here k = 27m /iV is a dimensionless wave vector and n is 
an integer. The wave vector is usually restricted to the 
interval —n<k<Tr but any other interval of length 2tt 
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can be chosen. The frequencies u> k are given by 



In the action and angle variables Hamiltonian has the 
form 
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Coefficient Vefc^fcV'fc'V'fc'" is given by 
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Here A/ is given in terms of Kronecker deltas as 

A/ = 2J ferro,! (6) 
m 

with to being an integer. Note that only terms with 
?7i = and to = ±1 have to be considered in (JSJ). Indeed, 
the maximum length for the sum of four wave vectors is 
4-7T but coefficient (j5jl vanishes in this case. 

In the problem of heat conduction the quantities of 
interest are the mode energy Ek = oJkJk and the total 
heat current given by 

jh = ^ V k UJ kJk = ^ V k E k, (?) 
k k 

where Vk is the group velocity. Note that Ek and jh 
represent only the harmonic parts of the corresponding 
quantities, it is assumed that the contributions from the 
anharmonic corrections are small for weak coupling. The 
approximate time evolution of the average energy of the 
normal mode for weak coupling and for the lattice with 
no temperature gradient and close to the thermal equi- 
librium is given by the homogeneous linearized Peierls 
equation 00,3- The equation is usually considered in 



the context of quantum mechanics for lattices with cubic 
anharmonicity. Derivation of this equation for the lattice 
with classical Hamiltonian (@J is straightforward. Note 
that in this case the following conditions on wave vectors 
and frequencies have to be satisfied simultaneously 

±k±k' ±k" ± k'" = 0, or ± 2tt, 

±LU ± J ± Lj" ± Uj" = (8) 

with the same ordering of signs for both relations. With 
the k dependence of frequencies © the relations can be 
satisfied only when two plus signs and two minus signs 
appear in ijHJl, i.e., in quantum mechanical terms, only the 
processes conserving the number of phonons contribute. 
In addition, although normal processes exist for this num- 
ber conserving case, they only result in exchange of quasi 
momenta between two colliding phonons, and, therefore, 
cannot change the phonon distribution. Thus, only umk- 
lapp number conserving processes contribute in the col- 
lision integral of the Peierls equation. We can write the 
linearized Peierls equation for the average energy Ek of 
mode fc, where overlining denotes averaging over a dis- 
tribution function. The equation has the form 
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The collision operator on the right hand side of equation inner product given by 
@ is a Hermitian operator in the Hilbert space with the 

(<7|/>= f dkglfk. (10) 

J — 7T 
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It can be shown that the collision operator has a con- 
tinuous spectrum that is bounded form below by zero 
[Ifil Il7j . We wrote equation (0 in the form that makes 
it easy to see that E k = const and E k = const /uj k are 
eigenstates of the collision operator with zero eigenval- 
ues. The first eigenstate corresponds to the conservation 
of the total energy. The second one corresponds to the 
conservation of the sum of action variables for all modes 
(or, in quantum mechanical language, to the conserva- 
tion of the number of phonons). Note that the second 
eigenstate has an infinite norm. 

We can write the average energy as 

E k = k B T + 6E k , (11) 

where k B T is the equilibrium value of E k and SE k is a 
deviation from that value. If we want the average ener- 
gies to approach their equilibrium value of fc^T for long 
times then 6E k should be orthogonal to both of the zero 
eigenvalue eigenstates of the collision operator Q, i.e., 
we must have 

dk5E k = 0, / dk — -=0. (12) 

-7T J-7T 
I 



The delta function appearing in JJjJ is meaningful only 
in the limit of N — > oo. In this limit we replace the 
sums by integrals and Kronecker deltas by delta function 
according to 



?f£^/*' ^h,k'^S(k-k>). (13) 



After the limit is taken we will have terms containing 
products of two delta functions in the integrand. There- 
fore, two integration can be rather easily performed. In- 
tegrations can be done much easier and the resulting ex- 
pressions have a simpler form if the wave vector is re- 
stricted to the interval from and 27r rather then from 
— 7r to 7r. Using explicit expressions for V kk i - k m and 
u) k we obtain after tedious but straightforward calcula- 
tions 
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Subscript int in the first integral in l|14|) means that the 
integral is taken over the interval where the integrand is 
real. This interval consists of two segments: from to 
l\{k) and from feik) to 2ir, where Zi(fc) and ^(fc) are the 
two solutions of the transcendental equation for k 1 



1 / k k'y k k' 

— I cos — + cos — I -sm-sm— = 0. (15) 



The solutions, which depend on A: as a parameter, satisfy 
h{k) < hik). In the next section we will use equation 
f|14|) to estimate the long time behavior of the heat cur- 
rent correlation function and the N dependence of the 
coefficient of thermal conductivity. 



III. THE LONG TIME BEHAVIOR OF THE 
CORRELATION FUNCTION 



relation function is defined as 



D N (t) = i j d{J k }d{a k }j{t)m e -^- 



(16) 



Here {J k } and {a k } denote the set of action and angle 
variables for all the modes, Z is the partition function 
for the equilibrium ensemble and j is the total energy 
current. The coefficient of heat conductivity is given by 

1 f* 

Kb-L t-*co J N-*oo 

We can rewrite D N (t) as 



D N {t) = - / d{J k }d{a k }$j(0)p(t) (18) 



with 



The coefficient of thermal conductivity can be calcu- 
lated by using the current correlation function. The cor- 
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where L is the Liouville operator corresponding to Hamil- 
tonian Q and £ is an auxiliary parameter insuring the 
correct dimensions for p(t). The parameter does not 
appear in the final expressions. In going from i|17|) 
to H8fl we also used the fact that the average current 
over the equilibrium distribution is zero. Equation l|18jl 
shows that the correlation function can be expressed 
through the average current per particle with the av- 
eraging performed over the nonequilibrium distribution 
function l|19l) . If we approximate the total heat current 
by its harmonic part J7J we can see that it depends only 
on the action variables and, therefore, the time evolution 
of the correlation function can be reduced to the time 
evolution of the average mode energies which is governed 
by Ijl4(l . Note that if only harmonic terms are kept in dis- 
tribution function l|19|) at t = then the initial average 
energy for mode k is 



E k (0) = k B T- 



2v k k 2 B T 2 



(20) 



This has the form given in Ijlljl with the deviation from 
ksT orthogonal to both of the zero eigenvalue eigenstates 
of the collision operator. Therefore, we expect the aver- 
age mode energies to approach fcsT for long times. To 
estimate the time behavior of E k based on l|14(l we will 
use the relaxation time approximation |3j. We assume 
that the energy of each mode approaches zero with a 
characteristic time T k which depends on the wave vector, 



i.e., 
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Some plausibility arguments in support of this approxi- 
mation were given in 0, Hfll | . In this approximation and 
for N — > oo the correlation function (for which we now 
drop the subscript N) is given by 
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Since the decay rate for energy of the normal mode with 
k = is zero (due to the conservation of the total mo- 
mentum) we can expect 1 /r^ to behave as some positive 
power of k for small k. The long time behavior of D(t) in 
()22J) will be determined by the small k behavior of 1/rfc. 
Following reference [f| we will further assume that the 
k dependence of l/r k for small k is the same as in the 
multiplicative part of the collision operator in (|14H . 

Both the relaxation time approximation and the as- 
sumption that the k dependence of the relaxation rate 
for small k is the same as in the multiplicative part of 
the collision operator has been widely used in the theory 
of heat conduction in insulators 0,0- A convincing justi- 
fication of both assumptions , however, is lacking. Refer- 
ence tries to justify both assumptions at least for wave 
vectors with small k by the following reasoning. If only 
the multiplicative part was kept in the collision operator 



the resulting equation would describe a physical situa- 
tion when all modes except mode k are in equilibrium. 
In general this is not the case. However, for any initial 
nonequilibrium distribution of energy all modes except 
those with very small k quickly relax to equilibrium. As 
a result as far as the small k modes are concerned after a 
short time the physical situation is similar to the one just 
described and the integral part of the collision operator 
becomes negligible compared to the multiplicative part. 

If we accept both approximation for equation l|14|) then 
we expect l/r^ cx sin 2 (fc/2)I(fc) with 
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This integral can be reduced to an elliptic integral of the 
first kind through the substitution x = tan(fc'/4), 



where 
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Integral (|24(l can be reduced to the Legendre normal form 
and its k dependence can be expressed in terms of the k 
dependence for the roots of the forth order polynomial 
(|2*5|l 2(|. Since the calculations are rather involved and 
we are interested only in the small k behavior of I(k) we 
give here a less rigorous but simpler estimation that gives 
the same result for small k. We expand the coefficients in 
the polynomial in powers of k and keep the lowest order 
terms in front of each monomial to get 



I(k) cx 



dx 



|| x 4 + Ak x 3 



2 X 
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Note that for positive k the denominator remains positive 
in the integration range since 



fr 4 k 2 
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(27) 



as can be checked by solving the corresponding quadratic 
equation for x 2 . Introducing the new variable y — k x ^x 
we obtain 



I(k) cx 



1 



/>oo 

x / dy- 

/0 Jh^lyi+^yS- + 4fc 2/3 y + 2 



(28) 



5 



The integral appearing in (|28|l is finite and remains finite 
for fc — 0. As a result, for small k we have /(fc) oc fc -1 / 3 
and, therefore, 

— oc fc 5 / 3 . (29) 

With this fc dependence for the relaxation rate the time 
dependence of the correlation function is determined by 
the following integral 

D(t)<x f dke- k5/3Kt vl (30) 
Jo 

Here if is a positive constant. Keeping in mind that Vk 
is a constant for small fc the long time behavior of D{t) 
is estimated to be 

D(t) oc (31) 

This implies that the heat conductivity coefficient n di- 
verges. Indeed, we have 

k oc lim / dr—— oc lim t 2/5 . (32) 

Clearly the divergence of k does not mean that the en- 
ergy propagates through the lattice instantaneously. It 
just implies that the Fourier heat law is not valid in the 
infinite FPU (3 lattice. We can also estimate the depen- 
dence of k on the size of the lattice. Following reference 
Q we restrict the integral in l|17|l to times smaller than 
the characteristic time for the energy propagation N/vk- 
This leads to the following N dependence for k, 

k oc N 2/5 . (33) 

Thus we can see that k diverges in the thermodynamic 
limit which is consistent with l|32(l . 

Apart from the fc dependence of the decay rate for 
the mode energy, its temperature dependence is also of 
interest. It follows from equation Ijl4|l that as a function 
of temperature the decay rate is proportional to T 2 . 

The estimations we have obtained are in rather good 
agreement with the results of computer simulations. The 
divergence of the heat conductivity coefficient as N 11 with 
[i « 0.37 was observed in the numerical studies of the 
FPU P model H H |H E3- This result is very close 
to estimation (j33J). Similarly, the fc 5 / 3 dependence of 
the decay rate for mode energy was observed 9]. The 
temperature dependence of the decay rate reported in Q 
is very close to T 2 for weak coupling. 

In reference the same types of t dependence for 
the correlation function and N dependence for k were 
obtained by using the mode coupling theory. As a func- 
tion of temperature the decay rate obtained through the 
mode coupling theory behaves as T 4 / 3 in the limit of 
weak coupling and as T 1 / 4 for strong coupling 0- Thus, 
the temperature dependence in the weak coupling limit 
is different from our T 2 estimation. Note, however, that 



according to Ref. [9j the mode coupling results should be 
valid for strong coupling and on very long time scales. In 
general, the mode coupling theory as used in allows 
to make some general statements about the long time 
behavior of the current correlation function for a class of 
one-dimensional lattices while equation (|14|l gives a more 
detailed picture of the energy equipartition between the 
normal modes for the special case of the j3 FPU lattice 
for the weak coupling case. If solved numerically, equa- 
tion 1141) will allow for the quantitative comparison of 
the energy equipartition given by the Peierls equation to 
the one observed in computer simulations. We will not 
attempt here to analyze the relation between our result 
and the mode coupling theory although this point clearly 
deserves attention. 

In a recent publication |21| it is claimed that k should 
diverge with system size L as L 1 ^ 3 for all momentum con- 
serving one dimensional systems. So far the most careful 
computer simulation [TTj | fail to confirm this claim. At 
present, therefore, this issue remains unsettled. 

It is well known that for systems such as a gas of hard 
spheres or Lorentz gas it is impossible to obtain the cor- 
rect long time behavior for the correlation functions if 
one uses only the kinetic equation p|. This is because 
for those systems the spectrum of the collision operator 
is discrete. As a result if only the kinetic equation is 
used the long time behavior is determined by the small- 
est non zero eigenvalue of the collision operator and has 
an exponential form. In contrast, in our case the collision 
operator has a continuous spectrum that is bounded form 
below by zero. This fact allows for the non trivial time 
dependence of the correlation function to be obtained al- 
ready in the framework of the kinetic equation. 



IV. CONCLUDING REMARKS 

Applying the Peierls equation to the FPU (3 lattice we 
estimated the wave vector and temperature dependence 
for the decay rate of the average mode energy, the long- 
time behavior of the current correlation function and the 
dependence of the coefficient of heat conductivity on the 
size of the lattice. The obtained results are in good agree- 
ment with the results of the recent computer simulations. 
As we used a number of strong assumptions it can be of 
interest to solve equation (|14fl numerically in order to ver- 
ify if the assumptions are correct and whether the time 
evolution of mode energies given by i|14|) is compatible 
with the results of computer simulations for the case of 
weak coupling. 

Recently lattices with external substrate potentials 
drew considerable attention since some of them show fi- 
nite heat conductivity for N —> oo |To L I2 I I2 I I2I |25 | . 
We can apply our analysis of Sec. Ill to show that the 
FPU lattice with added harmonic on-site potential of the 
form u 2 is likely to have finite heat conductivity for 
infinite lattice. It is easy to show that in this case for 
fc — > the harmonic frequency tends to a constant value 
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while the group velocity becomes proportional to k. The 
energy of the normal mode with k = is still a constant 
of motion since coefficient {5J vanishes when at least one 
the A;'s is zero. Therefore, we can expect the decay rate 
of the mode energy behave as k v for small k. This will 
lead to the i -3 / 1 ' long time behavior of the current cor- 
relation function and, therefore, finite heat conductivity 
for v < 3. Thus, if adding the harmonic on-site potential 
does not appreciably changes the /c 5//3 wave vector de- 
pendence of the decay rate, we can expect to find finite 
heat conductivity in this case. 
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